
Draft version August 29, 2005 

Preprint typeset using IAT^X style emulateapj v. 26/01/00 


EVIDENCE FOR HARMONIC CONTENT AND FREQUENCY EVOLUTION OF OSCILLATIONS 
DURING THE RISING PHASE OF X-RAY BURSTS FROM 4U 1636-536 

Sudip Bhattacharyya 1,2 , and Tod E. Strohmayer 2 
Draft version August 29, 2005 

ABSTRACT 

We report on a study of the evolution of burst oscillation properties during the rising phase of X-ray 
bursts from 4U 1636-536 observed with the proportional counter array (PCA) on board the Rossi X-Ray 
Timing Explorer (RXTE). We present evidence for significant harmonic structure of burst oscillation 
pulses during the early rising phases of bursts. This is the first such detection in burst rise oscillations, 
and is very important for constraining neutron star structure parameters and the equation of state 
models of matter at the core of a neutron star. The detection of harmonic content only during the initial 
portions of the burst rise is consistent with the theoretical expectation that with time the thermonuclear 
burning region becomes larger, and hence the fundamental and harmonic amplitudes both diminish. We 
also find, for the first time from this source, strong evidence of oscillation frequency increase during the 
burst rise. The timing behavior of harmonic content, amplitude, and frequency of burst rise oscillations 
may be important in understanding the spreading of thermonuclear flames under the extreme physical 
conditions on neutron star surfaces. 

Subject headings: equation of state — methods: data analysis — stars: neutron — X-rays: binaries — 
X-rays: bursts — X-rays: individual (4U 1636-536) 


1. INTRODUCTION 

Millisecond period brightness oscillations, “burst oscil- 
lations,” during thermonuclear (type I) X-ray bursts from 
the surfaces of neutron stars in low mass X-ray binary 
(LMXB) systems result from an asymmetric brightness 
pattern on the rotating stellar surface (Chakrabarty et al. 
2003; Strohmayer, & Bildsten 2003). This timing feature 
reveals the stellar spin period and can provide important 
information about the other stellar parameters (radius, 
mass, etc.; Miller, Sz Lamb 1998; Nath, Strohmayer, & 
Swank 2002; Muno, Ozel, & Chakrabarty 2002). Theo- 
retical modelling of these oscillations has the potential to 
constrain the equation of state (EOS) of the dense mat- 
ter at the core of a neutron star (Bhattacharyya et al. 
2005). This can be most effectively done if the burst oscil- 
lation has some harmonic content, as the fitting of a pure 
sinusoidal lightcurve does not put strong constraints on 
the stellar parameters. Although burst oscillations have 
been discovered from more than a dozen LMXBs, only 
one source (the accreting millisecond pulsar XTE J1814- 
338) has shown significant harmonic content (during burst 
decay; see Strohmayer et al. 2003; Watts, Strohmayer, 
& Markwardt 2005). In this Letter we focus on oscilla- 
tions during burst rise, because at the beginning of the 
burst, the size of the burning region (hot spot) is theo- 
retically expected to be at its smallest, and the harmonic 
content should theoretically be larger, and hence might be 
detected. 

The study of burst oscillations during burst rise is im- 
portant for another reason. At the onset of the burst, a 
small portion of the stellar surface is ignited, and then 
the burning region spreads and engulfs the whole stel- 
lar surface during the burst rise (Spitkovsky, Levin, & 


Ushomirsky 2002). This is natural, because simultane- 
ous ignition of the whole stellar surface would require very 
fine tuning. Understanding this spreading is important, as 
it involves complex nuclear physics and geophysical fluid 
dynamics (with significant stellar spin), under conditions 
of extreme gravity, magnetic field and radiation pressure. 
Other than the burst rise lightcurve, the time evolution of 
three properties (frequency, amplitude, and harmonic con- 
tent) of burst oscillations are the most useful observational 
tools to investigate this problem. Studies to date have not 
found harmonic content during burst rise. Strohmayer, 
Zhang, & Swank (1997) found evidence for a decrease in 
amplitude during burst rise from an analysis of RXTE data 
from the LMXB 4U 1728-34. Frequency increase during 
the burst rise from the accreting millisecond pulsar SAX 
J1808.4-3658 has also been reported (Chakrabarty et al. 
2003). 

In this Letter, we study burst oscillation properties of 
the LMXB 4U 1636-536 during burst rise. This is a regu- 
lar bursting source and the Rossi X-ray Timing Explorer 
(RXTE) has discovered oscillations during many bursts 
from this source over the course of its mission. The oscil- 
lation frequency is ~ 582 Hz. Here, we report the finding 
of frequency evolution during the rise of several bursts, as 
well as evidence for significant harmonic content of burst 
oscillations at the beginning of burst rise. 

2. DATA ANALYSIS AND RESULTS 

More than a hundred bursts have so far been observed 
from the LMXB 4U 1636-536 by the proportional counter 
array (PCA) on board RXTE. By analysing the archival 
data we find that 23 of them show at least moderately 
strong oscillations during burst rise. In order to search 
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sensitively for harmonic content, it is necessary to maxi- 
mize the fundamental power by taking any frequency evo- 
lution into account (see Miller 1999; Strohmayer & Mark- 
war dt 1999; Muno et al. 2002). Therefore, first we ex- 
plore frequency evolution during burst rise using the fol- 
lowing procedures: (1) we calculate dynamic power spec- 
tra (Strohmayer & Markwardt 1999) using a time duration 
that is small enough to resolve the burst rise, but is still 
large enough to accumulate significant signal power. The 
resulting dynamic spectra (see for example, panel a of Fig. 
1, and Fig. 2) provide an indication of the frequency evo- 
lution behavior. (2) We carry out a phase timing analysis 
(Muno et al. 2000) to confirm the indications in the dy- 
namic spectra. We divide the burst rise time interval into 
several bins of a fixed chosen length, and then assuming a 
frequency evolution model, we calculate the average phase 
(i/ik) in each bin ( k ). We use simple polynomial models to 
describe the frequency evolution. The corresponding x 2 is 
calculated using the formula \ 2 = “ ^k) 2 l G % k 

(Strohmayer h, Markwardt 2002), where M is the num- 
ber of bins, and a^ k = l/y/(Z 2 ) (Z\ is the fundamental 
power in the corresponding bin). We find the best fit pa- 
rameter values for a particular frequency evolution model 
by minimizing this \ 2 a nd we calculate the uncertainty in 
each parameter by increasing the \ 2 value by the appro- 
priate amount (Strohmayer Sz Markwardt 2002; Press et 
al. 1992). (3) To confirm these results, we calculate the 
total fundamental power ( Z \ ) during the burst rise time 
interval using these best fit frequency evolution model pa- 
rameter values, and ensure that this power is close to the 
maximum power obtained from any parameter values. 

We discover four bursts with significant frequency evo- 
lution during the rise. These are listed in Table 1. Burst 
A shows the clearest evidence of frequency increase (see 
panel a of Fig. 1). For this burst, a constant frequency 
model gives a reduced x 2 — 195.03/6 and fundamental 
power Z\ = 57.41. We also calculated the phase residuals 
(see Strohmayer &; Markwardt 2002 for details) for this 
model. The large systematic deviations from the mean 
value (panel b of Fig. 1), and the very high x 2 indicate that 
a constant frequency model for burst A can be strongly re- 
jected. Next, we add a linear term to the frequency model, 
and find that the corresponding reduced x 2 = 33.76/5, and 
Z\ = 134.15. This is a better fit, but still not statistically 
acceptable. Therefore, we include a nonlinear term to the 
frequency model. This model has a reduced x 2 = 2.00/4, 
and Z\ = 170.42, and the corresponding phase residuals 
have small random deviations from the mean value (panel 
c of Fig. 1). Hence, we conclude, that for burst A, the data 
strongly indicate a nonlinear frequency increase (Table 1). 

We fit constant frequency models to the other bursts 
listed in Table 1. For bursts B, C, & D, the reduced x 2 
(Zf) values for this model are 225.58/6 (86.91), 43.84/6 
(86.16), & 142.11/5 (20.42) respectively. These are all 
poor fits, and hence the constant frequency model for these 
bursts can be strongly rejected. Next, we include a linear 
term in the frequency models for these bursts. The cor- 
responding reduced x 2 (Z 2 ) values are 7.16/5 (200.35), 
12.70/5 (173.18), 8z 5.11/4 (89.75) respectively. These 
fits are acceptable for bursts B & D. For burst C, the 
high power value and the good visual fit of this frequency 
model to the power contours (panel b of Fig. 2) suggest 


that this model is on average correct, and the high reduced 
X 2 value may be caused by fluctuation of the frequency on 
short time scales. Therefore, for bursts B, C, &: D, a linear 
frequency increase (see Table 1) is highly probable. 

In panel d of Fig. 1, we show the rms amplitude vari- 
ation with time during the rise of burst A. The funda- 
mental amplitude shows an initially fast and then a slower 
decrease. Moreover, there is some weak indication of a 
significant 1st harmonic amplitude during the first ~ 1/2 
second of this burst. 

Using our best fitting frequency evolution models, we 
searched for harmonic content in the burst oscillations. 
Individually, none of the bursts shows strong harmonic 
power, therefore, we added the bursts coherently to get 
more signal. For this purpose, we chose the bursts with 
fundamental Z 2 power > 30 (for the constant frequency 
model) during burst rise, and added them (nine in num- 
ber; mentioned in Table 2) coherently (i.e., we shifted the 
phases suitably in order to maximize the total fundamen- 
tal power). The total co-added, phase-folded lightourve 
does not give a significant harmonic power. A possible 
reason for this may be that during a significant portion 
of the burst rise interval, the size of the burning region is 
large enough that the harmonic amplitude is small and not 
detectable. To further address this possibility we consider 
five time intervals (starting at the time of burst onset) 
of length; l/4th, l/3rd, half, 2/3rd, and all of the rise 
time. For each interval fraction, and for each burst, we 
consider that fraction of the total rise time, fit it with our 
frequency evolution model to get best fit parameter values, 
and use these best fit values to calculate the phases. We do 
this for all the nine bursts and then add them coherently 
(separately for each interval fraction). We find that the 
interval comprising 1 /3rd of the rise time gives the high- 
est harmonic power (17.52). The search at the harmonic 
frequency in any independent subinterval is essentially a 
single trial search (see Miller 1999), so we can estimate 
the significance of this value using a x 2 distribution with 2 
degrees of freedom. This gives a single-trial probability of 
1.57 x 10~ 4 to find a power as high by chance. We searched 
5 intervals, but they are not all independent, so the number 
of trials is between 1 and 5. Using 5 to get a conservative, 
lower bound on the probability gives 7.85 x 10 -4 , which 
is still a bit better than a 3cr detection. As we increase 
the time interval, (that is, use more of the rise), the har- 
monic power decreases gradually, which is consistent with 
the expectation that with time the burning region becomes 
larger, and hence the harmonic content diminishes. Due 
to the consistency with this theoretical expectation, and 
the better than 3 a harmonic power for l/3rd of the rise 
interval, we conclude that harmonic content in the burst 
rise oscillations has been marginally detected. 

We extracted the corresponding (i.e., for l/3rd of the 
rise time interval) phase-folded lightcurve from the data 
(after subtraction of the persistent emission) of the nine 
bursts. We fit it with two models: (1) a single sinusoid (the 
fundamental) around a constant level, and (2) two sinu- 
soids (fundamental and 1st harmonic) around a constant 
level. The former one gives a reduced x 2 value 26.26/13, 
while this value for the latter model is 8.33/11 (Table 2). 
These results support our finding that the addition of a 1st 
harmonic provides a better description of the data. Fig. 
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3 shows the data, the best fit model (model 2 of Table 2), 
and all the components of the model. 

3. DISCUSSION 

In this Letter, we report two new observational results: 
(1) the first detection of frequency evolution (increase) of 
burst rise oscillations from 4U 1636-536, and (2) the first 
evidence for harmonic content of burst rise oscillations 
(from any source). These effects may be a direct result 
of thermonuclear flame propagation on the slellar surface 
(Bhattacharyya & Strohmayer 2005a; 2005b). For exam- 
ple, consider ignition of a burst off of the equator in the 
northern hemisphere. Initially the burning region is rela- 
tively small. Its size is constrained by Coriolis forces and 
is given by the Rossby adjusment radius (Pedlosky 1987; 
Spitkovsky et al. 2002). This confined burning region can 
produce both a large fundamental and harmonic ampli- 
tude, more or less consistent with the observations. The 
frequency, when first observed, is at its lowest value and 
then increases monotonically. At least two effects can ac- 
count for this behavior; hydrostatic expansion — and sub- 
sequent spin- down — makes the burning region slip west- 
ward (on a star rotating eastward; see Strohmayer, Ja- 
hoda, Giles, & Lee 1997; Cumming & Bildsten 2001; hum- 
ming et al. 2002; Spitkovsky et al. 2002; Bhattacharyya 
& Strohmayer 2005b), and the southbound front will slip 
even further westward due to conservation of angular mo- 
mentum (Bhattacharyya & Strohmayer 2005b). Thus, the 
hot region initially has a net retrograde drift in the frame 
of the neutron star, so the observed frequency is less than 
the spin frequency. As the front approaches the equator 
it spreads faster, eventually forming a more or less sym- 
metric equatorial belt (Spitkovsky et al. 2002). This can 
plausibly reduce the pulsation amplitude in both the fun- 
damental and harmonic, and the westward drift slows be- 
cause mass elements moving southward below the equa- 
tor now drift eastward, conserving angular momentum. 
Thus, the frequency increases from its initial (low) value. 
Once the equatorial belt has been ignited, it seems likely 
that residual asymmetry associated with the initial, north- 
bound burning front is responsible for the observed, lower 
amplitude oscillations. Detailed calculations of the flux 
from such a spreading burning front are required to ex- 
plore this scenario quantitatively, but the discussion above 
provides a plausible qualitative understanding of many of 
the observed properties (more detailed discussions are in 
Bhattacharyya &; Strohmayer 2005b). 

As noted above, the harmonic content found at the be- 


ginning of the bursts is consistent with the expected small 
size of the burning region. Although in this study we do 
not fit the phase-folded lightcurve (solid curve in Fig. 3) 
with realistic theoretical model lightcurves (such fitting 
may not be so straightforward to interpret, as the burning 
region geometries for bursts may be different), we have 
computed theoretical models to show that the inferred 
harmonic content is theoretically possible. For example, 
the harmonic to fundamental amplitude ratio (a 2 /ai), and 
the relative phase difference (e x — e 2 ) of the components 
can be reproduced reasonably well with a model assuming 
emission from a circular hot spot (Bhattacharyya et al. 
2005). Using a dimensionless stellar radius-to-mass ratio 
Rc 2 /GM = 4.5, stellar mass M = I.SMq, observer’s incli- 
nation angle i = 70°, 0-position of the center of the circu- 
lar burning region 9 C — 50°, angular radius of the burning 
region A 9 — 25°, beaming parameter n = 1.0, and black- 
body temperature of the burning region Tbb = 2.0 keV, 
we can explain the relative strength and phasing of the 
fundamental and harmonic components. Here, we have 
assumed a Schwarzschild spacetime, and the beaming pa- 
rameter n gives a measure of the beaming in the frame 
corotating with the star (see Bhattacharyya et al. 2005 for 
a more detailed discussion of the model parameters). For 
the lightcurve calculated with these model parameter val- 
ues, we find a 2 /a x = 0.24 and e x — e 2 = 0.41, while for the 
data (Fig. 3), a 2 /a x = 0.24±0.06 and e x — e 2 = 0.36=b0.02. 
The constant level, required by the observed lightcurve, 
can plausibly be supplied by a symmetric belt-like compo- 
nent of the burning region (as mentioned earlier). 

Although we have found reasonably strong evidence for 
harmonic structure at the onset of bursts, it required the 
addition of a significant number of bursts to increase the 
signal to noise ratio. These results suggest that detailed 
studies of the onset of bursts can, in principle, provide 
insight into the structure of neutron stars, and how ther- 
monuclear flames propagate on their surfaces. Unfortu- 
nately, current studies are observationally limited by the 
detected count rates. For a larger count rate, the observed 
oscillation power in a given time interval would be greater. 
This would facilitate the detection of harmonic content, 
and also enable us to do more detailed studies of frequency 
and amplitude evolution. Therefore, new openings in this 
field are very likely with future larger area detectors. 


This work was supported in part by NASA Guest Inves- 
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Fig. 1. — Time evolution of different observed burst properties during the rise of burst A (see Table 1) from 4U 1636-536. Panel a gives 
the detected intensity (histogram), power contours (minimum and maximum power values are 18 and 55) using dynamic power spectra (for 
0.4 s duration at 0.02 s intervals), and the best fit model from Table 1. Panel b gives the phase residuals (phase varies from 0 to 1) and 
rms deviation (broken horizontal lines) for constant frequency (580.72 Hz) fitting. Panel c is same as panel 6, but for the best fit values of 
frequency evolution parameters (Table 1). Panel d gives the rms amplitudes of brightness oscillation (persistent emission subtracted). The 
upper histogram is for the fundamental amplitude, and the lower histogram is for the 1st harmonic amplitude. Here the horizontal lines give 
the binsize and the vertical lines give ler error. 
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Fig. 2. — Similar as panel a of Fig. 1 (panel a is for burst B, panel b is for burst C, and panel c is for burst D). For each of the panels, the 
dynamic power spectra are calculated for 0.3 s duration at 0.015 s intervals. The minimum and maximum power values of the contours are 
(23, 77), (20, 53), and (15, 46) for panels a , b and c respectively. 
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Phase 


Fig. 3. — Phase-folded lightcurve of burst oscillation during burst rise from 4U 1636-536. The solid curve shows the data (combined for 
nine bursts for l/3rd of the rise time interval; see Table 2), dashed curve is model 2 (see Table 2), solid horizontal line is the constant level 
of the model, dotted curve is the fundamental component of the model, and dash-dot curve is the 1st harmonic component of the model. 
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Table 1 

Frequency evolution model parameters*" (with la error) for four bursts. 


Obsld 

Start date 

Burst 

^0 

V 

V 

x 2 

Z\ h 

60032-05-03-00 

2002 Jan 12 

A 

577.70 

9 ^ q + 0.24 

-0.42l°;°7 

2.00 

170.42 

60032-05-03-00 

2002 Jan 13 

B 

579.19t°$ 

0 81 + 0-06 
U * Oi - 0.05 

- 

7.16 

200.35 

60032-05-05-00 

2002 Jan 14 

C 

580.62ioii2 

00 00 
o o 
d o 
+ 1 
00 
CN 
O 

- 

12.70 

173.18 

60032-05-10-00 

2002 Jan 22 

D 

578.44±°; 1 j 

1 21 +0 ' 11 

- 

5.11 

89.75 


a Frequency evolution model: v(t) = i/ 0 + vt + vt 1 , where i/ 0 — i/(0). 
b Fundamental power during burst rise. 
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Table 2 

Fitting of phase-folded lightcurve* : best fit parameter 15 values (with la error). 


Model 

ao 

a l 

Cl 


^2 

X 2 / do f 

i° 

667.87 ± 10.00 

241.42 ±14.23 

0.37 ±0.01 

- 

- 

26.26/13 

2 d 

668.99 ±10.00 

245.60 ±14.26 

0.37 ±0.01 

59.49 ±14.05 

0.01 ± 0.02 

8.33/11 


Combination (for first l/3rd of the rise time interval) of nine bursts from the following Obslds 
(start date): 30053-02-02-02 (1998 Aug 19), 30053-02-02-00 (1998 Aug 20), 40028-01-02-00 (1999 
Feb 27), 40028-01-08-00 (1999 Jun 18), 50030-02-01-00 (2000 Nov 05), 60032-05-03-00 (2002 Jan 
12), 60032-05-03-00 (2002 Jan 13), 60032-05-05-00 (2002 Jan 14) and 60032-05-13-00 (2002 Feb 
05). 

b Model lightcurve intensity is I = ao + ai sin(27r(e — ej)) + ct2 sin(47r(e — 62)), where e is the phase 
variable. 

c Only fundamental. 

d Fundamental + 1st harmonic. 



